Three Structure MOCHA

require(Ternary)

## Loading required package: Ternary

## Warning: package 'Ternary' was built under R version 4.3.2

Pcol <- rgb(123/255,188/255,73/255)
Vcol <- rgb(168/255,68/255,151/255)
Mcol <- rgb(231/255,186/255,96/255)

# Parameters

# Allocations
alphaV <- 0.3
alphaP <- 0.1
alphaM <- 1-alphaV-alphaP

# External resource availability
L <- 100
I <- 10
B <- 1e6

# Respiratory costs
r <- 0.1 #0.01

# Uptakes
yGM <- 1
uCL <- 0.01
uNI <- 1
uCB <- 0.000003
uNB <- uCB/10

# Stoichiometries
qC <- 1
qN <- 1
qG <- 1




# Set up the simulation
tset <- seq(from = 0, to = 10, by = 0.0001)
Vset <- NaN*tset; Vset[1] <- alphaV
Pset <- Vset; Pset[1] <- alphaP
Mset <- Vset; Mset[1] <- alphaM
growth <- NaN*tset

# Run the simulation
for(i in 2:length(tset)){
  dt <- tset[i] - tset[i-1]
  jC <- dt*( uCL*L*Pset[i-1] + uCB*B*Vset[i-1] - r*(Pset[i-1]+Vset[i-1]+Mset[i-1]) )
  jN <- dt*( uNI*I*Pset[i-1] + uNB*B*Vset[i-1] )
  jG <- dt*( yGM*Mset[i-1] )
  
  growth[i] <- min(qC*jC,qN*jN,qG*jG)
  
  Pset[i] <- Pset[i-1] + growth[i]*alphaP
  Vset[i] <- Vset[i-1] + growth[i]*alphaV
  Mset[i] <- Mset[i-1] + growth[i]*alphaM
    
    
    
}

plot(tset,Pset,type='l',las=1,xlab='Time',ylab='Organelle Structures',ylim=c(0,max(max(Pset),max(Mset),max(Vset))),col=Pcol,lwd=2)
lines(tset,Vset,lwd=2,col=Vcol)
lines(tset,Mset,lwd=2,col=Mcol)

plot(tset,alphaP*growth/Pset,type='l',las=1,col=Pcol,lwd=2,ylim=c(0,.15))
lines(tset,alphaM*growth/Mset,col=Mcol,lwd=2)
lines(tset,alphaV*growth/Vset,col=Vcol,lwd=2)

# Growth rates converge
lm(log(tail(Mset,50))~tail(tset,50))

## 
## Call:
## lm(formula = log(tail(Mset, 50)) ~ tail(tset, 50))
## 
## Coefficients:
##    (Intercept)  tail(tset, 50)  
##        -0.5108          0.6000

lm(log(tail(Vset,50))~tail(tset,50))

## 
## Call:
## lm(formula = log(tail(Vset, 50)) ~ tail(tset, 50))
## 
## Coefficients:
##    (Intercept)  tail(tset, 50)  
##         -1.204           0.600

lm(log(tail(Pset,50))~tail(tset,50))

## 
## Call:
## lm(formula = log(tail(Pset, 50)) ~ tail(tset, 50))
## 
## Coefficients:
##    (Intercept)  tail(tset, 50)  
##         -2.303           0.600

# Ratios converge on allocations
tail(Pset/(Pset+Mset+Vset))

## [1] 0.1 0.1 0.1 0.1 0.1 0.1

tail(Mset/(Pset+Mset+Vset)) 

## [1] 0.6 0.6 0.6 0.6 0.6 0.6

# Can we predict growth rate?
jC <- uCL*L*alphaP + uCB*B*alphaV - r
jN <- uNI*I*alphaP + uNB*B*alphaV
jG <- yGM*alphaM

guess_growth <- min(qC*jC,qN*jN,qG*jG)
guess_growth

## [1] 0.6

# It works!!

Show an example of convergence.

# Parameters

# Allocations
alphaV <- 0.3
alphaP <- 0.1
alphaM <- 1-alphaV-alphaP

# External resource availability
L <- 100
I <- 10
B <- 1e6

# Respiratory costs
r <- 0.1#0.01

# Uptakes
yGM <- 1
uCL <- 0.01
uNI <- 1
uCB <- 0.000003
uNB <- uCB/10

# Stoichiometries
qC <- 1
qN <- 1
qG <- 1




# Set up the simulation
tset <- seq(from = 0, to = 10, by = 0.0001)
Vset <- NaN*tset; Vset[1] <- 8
Pset <- Vset; Pset[1] <- 3
Mset <- Vset; Mset[1] <- 3.5
growth <- NaN*tset

# Run the simulation
for(i in 2:length(tset)){
  dt <- tset[i] - tset[i-1]
  jC <- dt*( uCL*L*Pset[i-1] + uCB*B*Vset[i-1] - r*(Pset[i-1]+Vset[i-1]+Mset[i-1]) )
  jN <- dt*( uNI*I*Pset[i-1] + uNB*B*Vset[i-1] )
  jG <- dt*( yGM*Mset[i-1] )
  
  growth[i] <- min(qC*jC,qN*jN,qG*jG)
  
  Pset[i] <- Pset[i-1] + growth[i]*alphaP
  Vset[i] <- Vset[i-1] + growth[i]*alphaV
  Mset[i] <- Mset[i-1] + growth[i]*alphaM
    
    
    
}


Tset <- Pset+Vset+Mset

plot(tset,Pset/Tset,col=Pcol,lwd=2,ylim=c(0,1),type='l',las=1,xlab='',ylab='')
abline(h=alphaP,lty=2,col=Pcol)
lines(tset,Vset/Tset,col=Vcol,lwd=2)
abline(h=alphaV,col=Vcol,lty=2)
lines(tset,Mset/Tset,col=Mcol,lwd=2)
abline(h=alphaM,col=Mcol,lty=2)
mtext('Time',side=1,line=2,cex=0.9)
mtext('Proportion of Biomass',side=2,line=2.5,cex=0.9)
legend(x='top',inset=0.05,legend=c('Vacuoles','Plastids','Machinery'),lwd=2,col=c(Vcol,Pcol,Mcol),horiz=T,cex=.8)

mtext(c('α_V','α_P','α_M'),side=4, line = 0.2, at = c(alphaV,alphaP,alphaM),las=1,cex=0.8)

Make a growth array

alphaPset <- seq(from=0, to = 1,by = 0.01)
alphaVset <- seq(from = 0, to = 1, by = 0.005)

growths <- matrix(data=NA,nrow=length(alphaPset),ncol=length(alphaVset))

for(i in 1:length(alphaPset)){
  for(j in 1:length(alphaVset)){
    if(alphaPset[i] + alphaVset[j] <=1){
      alphaP <- alphaPset[i]
      alphaV <- alphaVset[j]
      alphaM <- 1-alphaP-alphaV
      jC <- uCL*L*alphaP + uCB*B*alphaV - r
      jN <- uNI*I*alphaP + uNB*B*alphaV
      jG <- yGM*alphaM

      growths[i,j] <- min(qC*jC,qN*jN,qG*jG)
    }
  }
}
  

image(growths,xlab='alpha P', ylab='alpha V',las=1)

Make a few showing different outcomes.

alphaPset <- seq(from=0, to = 1,by = 0.01)
alphaVset <- seq(from = 0, to = 1, by = 0.05)

par(mar=c(4,4,1,1),mfrow=c(1,3))


# Option 1: Phototrophy
# External resource availability
L <- 100
I <- 5
B <- 0


growths <- matrix(data=NA,nrow=length(alphaPset),ncol=length(alphaVset))

for(i in 1:length(alphaPset)){
  for(j in 1:length(alphaVset)){
    if(alphaPset[i] + alphaVset[j] <=1){
      alphaP <- alphaPset[i]
      alphaV <- alphaVset[j]
      alphaM <- 1-alphaP-alphaV
      jC <- uCL*L*alphaP + uCB*B*alphaV - r
      jN <- uNI*I*alphaP + uNB*B*alphaV
      jG <- yGM*alphaM
      growths[i,j] <- min(qC*jC,qN*jN,qG*jG)
    }
  }
}

image(growths,xlab='alpha P', ylab='alpha V',las=1)



# Option 2: Phagotrophy
# External resource availability
L <- 0
I <- 0
B <- 1e6


growths <- matrix(data=NA,nrow=length(alphaPset),ncol=length(alphaVset))

for(i in 1:length(alphaPset)){
  for(j in 1:length(alphaVset)){
    if(alphaPset[i] + alphaVset[j] <=1){
      alphaP <- alphaPset[i]
      alphaV <- alphaVset[j]
      alphaM <- 1-alphaP-alphaV
      jC <- uCL*L*alphaP + uCB*B*alphaV - r
      jN <- uNI*I*alphaP + uNB*B*alphaV
      jG <- yGM*alphaM
      growths[i,j] <- min(qC*jC,qN*jN,qG*jG)
    }
  }
}

image(growths,xlab='alpha P', ylab='alpha V',las=1)



# Option 3: Mixotrophy
# External resource availability
L <- 100
I <- 5
B <- 1e6


growths <- matrix(data=NA,nrow=length(alphaPset),ncol=length(alphaVset))

for(i in 1:length(alphaPset)){
  for(j in 1:length(alphaVset)){
    if(alphaPset[i] + alphaVset[j] <=1){
      alphaP <- alphaPset[i]
      alphaV <- alphaVset[j]
      alphaM <- 1-alphaP-alphaV
      jC <- uCL*L*alphaP + uCB*B*alphaV - r
      jN <- uNI*I*alphaP + uNB*B*alphaV
      jG <- yGM*alphaM
      growths[i,j] <- min(qC*jC,qN*jN,qG*jG)
    }
  }
}

#image(growths,xlab='alpha P', ylab='alpha V',las=1, col=hcl.colors(30,'YlOrRd',rev=T))
image(growths,xlab='alpha P', ylab='alpha V',las=1)

OK so looking at just the mixotroph plot we can see that there are three “ridges” to the surface.

par(mar=c(4,4,1,1))

alphaVset <- seq(from = 0, to = 1, length.out = 500)
alphaPset <- seq(from = 0, to = 1, length.out = 495)

# Option 3: Mixotrophy
# External resource availability
L <- 100
I <- 5
B <- 1e6


growths <- matrix(data=NA,nrow=length(alphaPset),ncol=length(alphaVset))

for(i in 1:length(alphaPset)){
  for(j in 1:length(alphaVset)){
    if(alphaPset[i] + alphaVset[j] <=1){
      alphaP <- alphaPset[i]
      alphaV <- alphaVset[j]
      alphaM <- 1-alphaP-alphaV
      jC <- uCL*L*alphaP + uCB*B*alphaV - r
      jN <- uNI*I*alphaP + uNB*B*alphaV
      jG <- yGM*alphaM
      growths[i,j] <- min(qC*jC,qN*jN,qG*jG)
    }
  }
}

image(growths,xlab='alpha P', ylab='alpha V',las=1, col=hcl.colors(100,'YlOrRd',rev=T))

Let’s see if we can trace these.

par(mar=c(4,4,1,1),mfrow=c(1,2))
image(growths,xlab='alpha P', ylab='alpha V',las=1, col=hcl.colors(100,'YlOrRd',rev=T))


alphaPset.1 <- (qN*uNB*B*alphaVset-qC*(uCB*B*alphaVset-r))/(qC*uCL*L-qN*uNI*I)
alphaPset.2 <- (qG*yGM*(1-alphaVset)+qC*r-qC*uCB*B*alphaVset)/(qG*yGM+qC*uCL*L)
alphaPset.3 <- (qG*yGM*(1-alphaVset) - qN*uNB*B*alphaVset)/(qN*uNI*I + qG*yGM)


image(growths,xlab='alpha P', ylab='alpha V',las=1, col=hcl.colors(100,'YlOrRd',rev=T))
lines(alphaPset.2,alphaVset,lty=2)
lines(alphaPset.1,alphaVset)
lines(alphaPset.3,alphaVset,lty=3)

legend(x='topright',inset=0.02,legend=c('jC=jN','jC=jG','jN=jG'),lty=c(1,2,3))

Finding the intersection…

plot(alphaVset,abs(alphaPset.1-alphaPset.2),type='l',xlab='alpha V',ylab='delta alpha P')
lines(alphaVset,abs(alphaPset.2-alphaPset.3),col='dodgerblue')
lines(alphaVset,abs(alphaPset.1-alphaPset.3),col='darkgreen')

points(alphaVset[which(abs(alphaPset.1-alphaPset.2)==min(abs(alphaPset.1-alphaPset.2)))],0)

#alphaPset.1[which(abs(alphaPset.1-alphaPset.2)==min(abs(alphaPset.1-alphaPset.2)))]

L <- 10
I <- 0
B <- 1e6

growths <- matrix(data=NA,nrow=length(alphaPset),ncol=length(alphaVset))

for(i in 1:length(alphaPset)){
  for(j in 1:length(alphaVset)){
    if(alphaPset[i] + alphaVset[j] <=1){
      alphaP <- alphaPset[i]
      alphaV <- alphaVset[j]
      alphaM <- 1-alphaP-alphaV
      jC <- uCL*L*alphaP + uCB*B*alphaV - r
      jN <- uNI*I*alphaP + uNB*B*alphaV
      jG <- yGM*alphaM
      growths[i,j] <- min(qC*jC,qN*jN,qG*jG)
    }
  }
}


par(mar=c(4,4,1,1),mfrow=c(1,2))
image(growths,xlab='alpha P', ylab='alpha V',las=1, col=hcl.colors(100,'YlOrRd',rev=T))

alphaVset2 <- seq(from = 0, to = 1, length.out = 500)


alphaPset.1 <- (qN*uNB*B*alphaVset2-qC*(uCB*B*alphaVset2-r))/(qC*uCL*L-qN*uNI*I)
alphaPset.2 <- (qG*yGM*(1-alphaVset2)+qC*r-qC*uCB*B*alphaVset2)/(qG*yGM+qC*uCL*L)
alphaPset.3 <- (qG*yGM*(1-alphaVset2) - qN*uNB*B*alphaVset2)/(qN*uNI*I + qG*yGM)


image(growths,xlab='alpha P', ylab='alpha V',las=1, col=hcl.colors(100,'YlOrRd',rev=T))
lines(alphaPset.2,alphaVset2,lty=2)
lines(alphaPset.1,alphaVset2)
lines(alphaPset.3,alphaVset2,lty=3)

points(alphaPset.1[which(abs(alphaPset.1-alphaPset.2)==min(abs(alphaPset.1-alphaPset.2)))],alphaVset2[which(abs(alphaPset.1-alphaPset.2)==min(abs(alphaPset.1-alphaPset.2)))],pch=21,lwd=1,cex=1.5)

term.a <- (qN*uNB*B-uCB*qC*B)/(qC*uCL*L-qN*uNI*I)
term.b <- (qG*yGM+uCB*qC*B)/(qG*yGM+qC*uCL*L)
term.c <- (qG*yGM+qC*r)/(qG*yGM+qC*uCL*L)
term.d <- (qC*r)/(qC*uCL*L-qN*uNI*I)

alphaVstar <- (term.c-term.d)/(term.a+term.b)
alphaPstar <- (qG*yGM*(1-alphaVstar)+qC*r-qC*uCB*B*alphaVstar)/(qG*yGM+qC*uCL*L)

points(alphaPstar,alphaVstar,pch=21,col='red',cex=2.25,lwd=1)


legend(x='topright',inset=0.02,legend=c('gC=gN','gC=gG','gN=gG','g* numer.','g* analyt.'),lty=c(1,2,3,0,0),pch=1,pt.cex=c(0,0,0,1.5,2.25),col=c('black','black','black','black','red'))

alphaVstar

## [1] 0

alphaPstar

## [1] 1

alphaMstar <- 1 - alphaVstar - alphaPstar
alphaMstar

## [1] 0

dim(growths)

## [1] 495 500

gstar <- max(growths,na.rm = T)
which(growths==gstar)

## [1] 190081

which(growths==max(growths,na.rm = T), arr.ind=T)

##      row col
## [1,]   1 385

alphaVset[which(growths==max(growths,na.rm = T), arr.ind=T)[2]]

## [1] 0.7695391

alphaPset[which(growths==max(growths,na.rm = T), arr.ind=T)[1]]

## [1] 0

(alphaVstar-alphaVset[which(growths==max(growths,na.rm = T), arr.ind=T)[2]])/(alphaVstar)*100

## [1] Inf

Ternary plotting

# Plot parameters
par(mar=rep(0.2, 4),mfrow=c(1,3))
colfunc <- colorRampPalette(c("black",'gray65',"white"))
resolutionsetting <- 100L


# MIXOTROPH PARAMETERS
# Respiratory costs
r <- 0.01

# Uptakes
yGM <- 1
uCL <- 0.01
uNI <- 1
uCB <- 0.000003
uNB <- uCB/10

# Stoichiometries
qC <- 1
qN <- 1
qG <- 1

# Growth rate functions
jCfcn <- function (a, b, c) { qC*(uCL*L*a + uCB*B*c - r) }
jNfcn <- function(a,b,c){   qN*(uNI*I*a + uNB*B*c) }
jGfcn <- function(a,b,c){   qG*yGM*b }


# PHOTOTROPH
# Make the plot
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Phototrophy Optimal',col='black')

# External resource availability
L <- 100
I <- 5
B <- 0

# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){   values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }

# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))
#ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))


# PHAGOTROPH
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Phagotrophy Optimal',col='black')

# External resource availability
L <- 0
I <- 0
B <- 1e6

# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){   values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }

# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))
#ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))


# MIXOTROPH
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Mixotrophy Optimal',col='black')

# External resource availability
L <- 10
I <- 2
B <- 1e6

# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){   values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }

# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))

#ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))

Now that we can make these plots, can we add lines to them?

par(mar=rep(0.2, 4))
resolutionsetting <- 100

# MIXOTROPH
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Mixotrophy Optimal',col='black')

# External resource availability
L <- 10
I <- 2
B <- 1e6

# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){   values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }


# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))
ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))

# Add analytical solution
term.a <- (qN*uNB*B-uCB*qC*B)/(qC*uCL*L-qN*uNI*I)
term.b <- (qG*yGM+uCB*qC*B)/(qG*yGM+qC*uCL*L)
term.c <- (qG*yGM+qC*r)/(qG*yGM+qC*uCL*L)
term.d <- (qC*r)/(qC*uCL*L-qN*uNI*I)

alphaVstar <- (term.c-term.d)/(term.a+term.b)
alphaPstar <- (qG*yGM*(1-alphaVstar)+qC*r-qC*uCB*B*alphaVstar)/(qG*yGM+qC*uCL*L)
alphaMstar <- 1 - alphaVstar - alphaPstar


# Zero-waste lines
alphaVset2 <- seq(from = 0, to = 1, length.out = 100)

alphaPset.1 <- (qN*uNB*B*alphaVset2-qC*(uCB*B*alphaVset2-r))/(qC*uCL*L-qN*uNI*I)
alphaPset.2 <- (qG*yGM*(1-alphaVset2)+qC*r-qC*uCB*B*alphaVset2)/(qG*yGM+qC*uCL*L)
alphaPset.3 <- (qG*yGM*(1-alphaVset2) - qN*uNB*B*alphaVset2)/(qN*uNI*I + qG*yGM)

points.to.plot <- as.data.frame(cbind(alphaPset.1,(1-alphaPset.1-alphaVset2),alphaVset2))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,1]>0 & points.to.plot[,2]>alphaMstar,]

TernaryPoints(points.to.plot[,1:3],type='l',lwd=0.5)

points.to.plot <- as.data.frame(cbind(alphaPset.2,(1-alphaPset.2-alphaVset2),alphaVset2))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,1]>0 & points.to.plot[,2]>0&points.to.plot[,2]<alphaMstar,]

TernaryPoints(points.to.plot[,1:3],type='l',lwd=0.5,lty=2)

points.to.plot <- as.data.frame(cbind(alphaPset.3,(1-alphaPset.3-alphaVset2),alphaVset2))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,1]>0 & points.to.plot[,2]>0&points.to.plot[,2]<alphaMstar,]

TernaryPoints(points.to.plot[,1:3],type='l',lwd=0.5,lty=3)



TernaryPoints(c(alphaPstar,alphaMstar,alphaVstar),pch=21,col='red',cex=1,lwd=1)

How about a strict strategy?

par(mar=rep(0.2, 4))
resolutionsetting <- 100

# MIXOTROPH
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Phagotrophy Optimal',col='black')

# External resource availability
L <- 0
I <- 0
B <- 1e6

# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){   values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }


# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))
ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))

# Add analytical solution
alphaVset2 <- seq(from = 0, to = 1, length.out = 100)

alphaVstar.1 <- (qG*yGM+r*qC)/(qC*uCB*B+qG*yGM)
alphaVstar.2 <- qG*yGM/(qN*uNB*B+qG*yGM)
if(alphaVstar.1 > alphaVstar.2){ 
  alphaVstar <- alphaVstar.1
  alphaMset <- qC*(uCB*B*alphaVset2 - r)/(qG*yGM)
} else { 
    alphaVstar <- alphaVstar.2
    alphaMset <- qN*(uNB*B*alphaVset2)/(qG*yGM)
    }

alphaPstar <- 0
alphaMstar <- 1 - alphaVstar - alphaPstar


# Zero-waste lines

alphaPset <- 1-alphaMset-alphaVset2

points.to.plot <- as.data.frame(cbind(alphaPset,alphaMset,alphaVset2))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,1]>=0,]

TernaryPoints(points.to.plot[,1:3],type='l',lwd=0.5)




TernaryPoints(c(1-alphaMstar-alphaVstar,alphaMstar,alphaVstar),pch=21,col='red',cex=1,lwd=1)

And now let’s try marking all three.

# Plot parameters
par(mar=rep(0.2, 4),mfrow=c(1,3))
colfunc <- colorRampPalette(c("black",'gray65',"white"))
resolutionsetting <- 100L


# MIXOTROPH PARAMETERS
# Respiratory costs
r <- 0.01

# Uptakes
yGM <- 1
uCL <- 0.01
uNI <- 1
uCB <- 0.000003
uNB <- uCB/10

# Stoichiometries
qC <- 1
qN <- 1
qG <- 1

# Growth rate functions
jCfcn <- function (a, b, c) { qC*(uCL*L*a + uCB*B*c - r) }
jNfcn <- function(a,b,c){   qN*(uNI*I*a + uNB*B*c) }
jGfcn <- function(a,b,c){   qG*yGM*b }


# PHOTOTROPH
# Make the plot
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Phototrophy Optimal',col='black')

# External resource availability
L <- 100
I <- 5
B <- 0

# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){   values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }

# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))
ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))

# Add analytical solution
alphaPset2 <- seq(from = 0, to = 1, length.out = 100)

alphaPstar.1 <- (qG*yGM+r*qC)/(qC*uCL*L+qG*yGM)
alphaPstar.2 <- qG*yGM/(qN*uNI*I+qG*yGM)
if(alphaPstar.1 > alphaPstar.2){ 
  alphaPstar <- alphaPstar.1
  alphaMset <- qC*(uCL*L*alphaPset2 - r)/(qG*yGM)
  linetypesetting <- 2
} else { 
    alphaPstar <- alphaPstar.2
    alphaMset <- qN*(uNI*I*alphaPset2)/(qG*yGM)
    linetypesetting <- 3
    }

alphaVstar <- 0
alphaMstar <- 1 - alphaVstar - alphaPstar


# Zero-waste lines

alphaVset <- 1-alphaMset-alphaPset2

points.to.plot <- as.data.frame(cbind(alphaPset2,alphaMset,alphaVset))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,3]>=0,]

TernaryPoints(points.to.plot[,1:3],type='l',lwd=2,lty = linetypesetting)

TernaryPoints(c(alphaPstar,alphaMstar,alphaVstar),cex=1,col='red',bg='red',pch=21)



# PHAGOTROPH
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Phagotrophy Optimal',col='black')

# External resource availability
L <- 0
I <- 0
B <- 1e6

# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){   values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }

# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))
ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))

# Add analytical solution
alphaVset2 <- seq(from = 0, to = 1, length.out = 100)

alphaVstar.1 <- (qG*yGM+r*qC)/(qC*uCB*B+qG*yGM)
alphaVstar.2 <- qG*yGM/(qN*uNB*B+qG*yGM)
if(alphaVstar.1 > alphaVstar.2){ 
  alphaVstar <- alphaVstar.1
  alphaMset <- qC*(uCB*B*alphaVset2 - r)/(qG*yGM)
  linetypesetting <- 2
} else { 
    alphaVstar <- alphaVstar.2
    alphaMset <- qN*(uNB*B*alphaVset2)/(qG*yGM)
    linetypesetting <- 3
    }

alphaPstar <- 0
alphaMstar <- 1 - alphaVstar - alphaPstar


# Zero-waste lines

alphaPset <- 1-alphaMset-alphaVset2

points.to.plot <- as.data.frame(cbind(alphaPset,alphaMset,alphaVset2))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,1]>=0,]

TernaryPoints(points.to.plot[,1:3],type='l',lwd=2,lty = linetypesetting)

TernaryPoints(c(alphaPstar,alphaMstar,alphaVstar),cex=1,col='red',bg='red',pch=21)


# MIXOTROPH
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Mixotrophy Optimal',col='black')

# External resource availability
L <- 10
I <- 2
B <- 1e6

# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){   values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }

# Add shading to the plot
# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))
ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))

# Add analytical solution
term.a <- (qN*uNB*B-uCB*qC*B)/(qC*uCL*L-qN*uNI*I)
term.b <- (qG*yGM+uCB*qC*B)/(qG*yGM+qC*uCL*L)
term.c <- (qG*yGM+qC*r)/(qG*yGM+qC*uCL*L)
term.d <- (qC*r)/(qC*uCL*L-qN*uNI*I)

alphaVstar <- (term.c-term.d)/(term.a+term.b)
alphaPstar <- (qG*yGM*(1-alphaVstar)+qC*r-qC*uCB*B*alphaVstar)/(qG*yGM+qC*uCL*L)
alphaMstar <- 1 - alphaVstar - alphaPstar


# Zero-waste lines
alphaVset2 <- seq(from = 0, to = 1, length.out = 100)

alphaPset.1 <- (qN*uNB*B*alphaVset2-qC*(uCB*B*alphaVset2-r))/(qC*uCL*L-qN*uNI*I)
alphaPset.2 <- (qG*yGM*(1-alphaVset2)+qC*r-qC*uCB*B*alphaVset2)/(qG*yGM+qC*uCL*L)
alphaPset.3 <- (qG*yGM*(1-alphaVset2) - qN*uNB*B*alphaVset2)/(qN*uNI*I + qG*yGM)

points.to.plot <- as.data.frame(cbind(alphaPset.1,(1-alphaPset.1-alphaVset2),alphaVset2))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,1]>0 & points.to.plot[,2]>alphaMstar,]

TernaryPoints(points.to.plot[,1:3],type='l',lwd=2)

points.to.plot <- as.data.frame(cbind(alphaPset.2,(1-alphaPset.2-alphaVset2),alphaVset2))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,1]>0 & points.to.plot[,2]>0&points.to.plot[,2]<alphaMstar,]

TernaryPoints(points.to.plot[,1:3],type='l',lwd=2,lty=2)

points.to.plot <- as.data.frame(cbind(alphaPset.3,(1-alphaPset.3-alphaVset2),alphaVset2))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,1]>0 & points.to.plot[,2]>0&points.to.plot[,2]<alphaMstar,]

TernaryPoints(points.to.plot[,1:3],type='l',lwd=2,lty=3)

TernaryPoints(c(alphaPstar,alphaMstar,alphaVstar),cex=1,col='red',bg='red',pch=21)

And all three in grayscale.

# Plot parameters
par(mar=rep(0.2, 4),mfrow=c(1,3))
colfunc <- colorRampPalette(c("black",'gray65',"white"))
resolutionsetting <- 100L


# MIXOTROPH PARAMETERS
# Respiratory costs
r <- 0.01

# Uptakes
yGM <- 1
uCL <- 0.01
uNI <- 1
uCB <- 0.000003
uNB <- uCB/10

# Stoichiometries
qC <- 1
qN <- 1
qG <- 1

# Growth rate functions
jCfcn <- function (a, b, c) { qC*(uCL*L*a + uCB*B*c - r) }
jNfcn <- function(a,b,c){   qN*(uNI*I*a + uNB*B*c) }
jGfcn <- function(a,b,c){   qG*yGM*b }


# PHOTOTROPH
# Make the plot
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Phototrophy Optimal',col='black')

# External resource availability
L <- 100
I <- 5
B <- 0

# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){   values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }

# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))
#ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))

aP1 <- (qG*yGM+r*qC)/(qC*uCL*L+qG*yGM)
aP2 <- (qG*yGM)/(qN*uNI*I+qG*yGM)
alphaPstar <- max(aP1,aP2)
alphaVstar <- 0
alphaMstar <- 1 - alphaPstar

TernaryPoints(c(alphaPstar,alphaMstar,alphaVstar),cex=1,col='red',bg='red',pch=21)



# PHAGOTROPH
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Phagotrophy Optimal',col='black')

# External resource availability
L <- 0
I <- 0
B <- 1e6

# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){   values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }

# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))
#ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))


aV1 <- (qG*yGM+r*qC)/(qC*uCB*B+qG*yGM)
aV2 <- (qG*yGM)/(qN*uNB*B+qG*yGM)
alphaVstar <- max(aV1,aV2)
alphaPstar <- 0
alphaMstar <- 1 - alphaVstar

TernaryPoints(c(alphaPstar,alphaMstar,alphaVstar),cex=1,col='red',bg='red',pch=21)


# MIXOTROPH
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Mixotrophy Optimal',col='black')

# External resource availability
L <- 10
I <- 2
B <- 1e6

# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){   values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }

# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))
#ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))

term.a <- (qN*uNB*B-uCB*qC*B)/(qC*uCL*L-qN*uNI*I)
term.b <- (qG*yGM+uCB*qC*B)/(qG*yGM+qC*uCL*L)
term.c <- (qG*yGM+qC*r)/(qG*yGM+qC*uCL*L)
term.d <- (qC*r)/(qC*uCL*L-qN*uNI*I)

alphaVstar <- (term.c-term.d)/(term.a+term.b)
alphaPstar <- (qG*yGM*(1-alphaVstar)+qC*r-qC*uCB*B*alphaVstar)/(qG*yGM+qC*uCL*L)
alphaMstar <- 1 - alphaVstar - alphaPstar

TernaryPoints(c(alphaPstar,alphaMstar,alphaVstar),cex=1,col='red',bg='red',pch=21)

# Plot parameters
par(mar=rep(0.2, 4),mfrow=c(1,3))
colfunc <- colorRampPalette(c("black",'gray65',"white"))
resolutionsetting <- 100L


# MIXOTROPH PARAMETERS
# Respiratory costs
r <- 0.01

# Uptakes
yGM <- 1
uCL <- 0.01
uNI <- 1
uCB <- 0.000003
uNB <- uCB/10

# Stoichiometries
qC <- 1
qN <- 1
qG <- 1

# Growth rate functions
jCfcn <- function (a, b, c) { qC*(uCL*L*a + uCB*B*c - r) }
jNfcn <- function(a,b,c){   qN*(uNI*I*a + uNB*B*c) }
jGfcn <- function(a,b,c){   qG*yGM*b }


# PHOTOTROPH
# Make the plot
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Phototrophy Optimal',col='black')

# External resource availability
L <- 100
I <- 5
B <- 0

# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){   values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }

# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))
#ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))

# Add analytical solution
alphaPset2 <- seq(from = 0, to = 1, length.out = 100)

alphaPstar.1 <- (qG*yGM+r*qC)/(qC*uCL*L+qG*yGM)
alphaPstar.2 <- qG*yGM/(qN*uNI*I+qG*yGM)
if(alphaPstar.1 > alphaPstar.2){ 
  alphaPstar <- alphaPstar.1
  alphaMset <- qC*(uCL*L*alphaPset2 - r)/(qG*yGM)
  linetypesetting <- 2
} else { 
    alphaPstar <- alphaPstar.2
    alphaMset <- qN*(uNI*I*alphaPset2)/(qG*yGM)
    linetypesetting <- 3
    }

alphaVstar <- 0
alphaMstar <- 1 - alphaVstar - alphaPstar


# Zero-waste lines

alphaVset <- 1-alphaMset-alphaPset2

points.to.plot <- as.data.frame(cbind(alphaPset2,alphaMset,alphaVset))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,3]>=0,]

TernaryPoints(points.to.plot[,1:3],type='l',lwd=1,lty = linetypesetting)

TernaryPoints(c(alphaPstar,alphaMstar,alphaVstar),cex=1,col='red',bg='red',pch=21)



# PHAGOTROPH
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Phagotrophy Optimal',col='black')

# External resource availability
L <- 0
I <- 0
B <- 1e6

# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){   values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }

# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))
#ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))

# Add analytical solution
alphaVset2 <- seq(from = 0, to = 1, length.out = 100)

alphaVstar.1 <- (qG*yGM+r*qC)/(qC*uCB*B+qG*yGM)
alphaVstar.2 <- qG*yGM/(qN*uNB*B+qG*yGM)
if(alphaVstar.1 > alphaVstar.2){ 
  alphaVstar <- alphaVstar.1
  alphaMset <- qC*(uCB*B*alphaVset2 - r)/(qG*yGM)
  linetypesetting <- 2
} else { 
    alphaVstar <- alphaVstar.2
    alphaMset <- qN*(uNB*B*alphaVset2)/(qG*yGM)
    linetypesetting <- 3
    }

alphaPstar <- 0
alphaMstar <- 1 - alphaVstar - alphaPstar


# Zero-waste lines

alphaPset <- 1-alphaMset-alphaVset2

points.to.plot <- as.data.frame(cbind(alphaPset,alphaMset,alphaVset2))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,1]>=0,]

TernaryPoints(points.to.plot[,1:3],type='l',lwd=1,lty = linetypesetting)

TernaryPoints(c(alphaPstar,alphaMstar,alphaVstar),cex=1,col='red',bg='red',pch=21)


# MIXOTROPH
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Mixotrophy Optimal',col='black')

# External resource availability
L <- 10
I <- 2
B <- 1e6

# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){   values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }

# Add shading to the plot
# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))
#ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))

# Add analytical solution
term.a <- (qN*uNB*B-uCB*qC*B)/(qC*uCL*L-qN*uNI*I)
term.b <- (qG*yGM+uCB*qC*B)/(qG*yGM+qC*uCL*L)
term.c <- (qG*yGM+qC*r)/(qG*yGM+qC*uCL*L)
term.d <- (qC*r)/(qC*uCL*L-qN*uNI*I)

alphaVstar <- (term.c-term.d)/(term.a+term.b)
alphaPstar <- (qG*yGM*(1-alphaVstar)+qC*r-qC*uCB*B*alphaVstar)/(qG*yGM+qC*uCL*L)
alphaMstar <- 1 - alphaVstar - alphaPstar


# Zero-waste lines
alphaVset2 <- seq(from = 0, to = 1, length.out = 100)

alphaPset.1 <- (qN*uNB*B*alphaVset2-qC*(uCB*B*alphaVset2-r))/(qC*uCL*L-qN*uNI*I)
alphaPset.2 <- (qG*yGM*(1-alphaVset2)+qC*r-qC*uCB*B*alphaVset2)/(qG*yGM+qC*uCL*L)
alphaPset.3 <- (qG*yGM*(1-alphaVset2) - qN*uNB*B*alphaVset2)/(qN*uNI*I + qG*yGM)

points.to.plot <- as.data.frame(cbind(alphaPset.1,(1-alphaPset.1-alphaVset2),alphaVset2))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,1]>0 & points.to.plot[,2]>alphaMstar,]

TernaryPoints(points.to.plot[,1:3],type='l',lwd=1)

points.to.plot <- as.data.frame(cbind(alphaPset.2,(1-alphaPset.2-alphaVset2),alphaVset2))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,1]>0 & points.to.plot[,2]>0&points.to.plot[,2]<alphaMstar,]

TernaryPoints(points.to.plot[,1:3],type='l',lwd=1,lty=2)

points.to.plot <- as.data.frame(cbind(alphaPset.3,(1-alphaPset.3-alphaVset2),alphaVset2))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,1]>0 & points.to.plot[,2]>0&points.to.plot[,2]<alphaMstar,]

TernaryPoints(points.to.plot[,1:3],type='l',lwd=1,lty=3)

TernaryPoints(c(alphaPstar,alphaMstar,alphaVstar),cex=1,col='red',bg='red',pch=21)

# Plot parameters
par(mar=rep(0.2, 4),mfrow=c(1,3))
colfunc <- colorRampPalette(c("black",'gray65',"white"))
resolutionsetting <- 100L


# MIXOTROPH PARAMETERS
# Respiratory costs
r <- 0.01

# Uptakes
yGM <- 1
uCL <- 0.01
uNI <- 1
uCB <- 0.000003
uNB <- uCB/10

# Stoichiometries
qC <- 1
qN <- 1
qG <- 1

# Growth rate functions
jCfcn <- function (a, b, c) { qC*(uCL*L*a + uCB*B*c - r) }
jNfcn <- function(a,b,c){   qN*(uNI*I*a + uNB*B*c) }
jGfcn <- function(a,b,c){   qG*yGM*b }


# PHOTOTROPH
# Make the plot
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Phototrophy Optimal',col='black')

# External resource availability
L <- 100
I <- 5
B <- 0

# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){   values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }

# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))
#ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))

# Add analytical solution
alphaPset2 <- seq(from = 0, to = 1, length.out = 100)

alphaPstar.1 <- (qG*yGM+r*qC)/(qC*uCL*L+qG*yGM)
alphaPstar.2 <- qG*yGM/(qN*uNI*I+qG*yGM)
if(alphaPstar.1 > alphaPstar.2){ 
  alphaPstar <- alphaPstar.1
  alphaMset <- qC*(uCL*L*alphaPset2 - r)/(qG*yGM)
  linetypesetting <- 2
} else { 
    alphaPstar <- alphaPstar.2
    alphaMset <- qN*(uNI*I*alphaPset2)/(qG*yGM)
    linetypesetting <- 3
    }

alphaVstar <- 0
alphaMstar <- 1 - alphaVstar - alphaPstar


# Zero-waste lines

alphaVset <- 1-alphaMset-alphaPset2

points.to.plot <- as.data.frame(cbind(alphaPset2,alphaMset,alphaVset))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,3]>=0,]

TernaryPoints(points.to.plot[,1:3],type='l',lwd=1,lty = linetypesetting)

TernaryPoints(c(alphaPstar,alphaMstar,alphaVstar),cex=1,col='red',bg='red',pch=21)



# PHAGOTROPH
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Phagotrophy Optimal',col='black')

# External resource availability
L <- 0
I <- 0
B <- 0.88e6 #1e6

# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){   values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }

# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))
#ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))

# Add analytical solution
alphaVset2 <- seq(from = 0, to = 1, length.out = 100)

alphaVstar.1 <- (qG*yGM+r*qC)/(qC*uCB*B+qG*yGM)
alphaVstar.2 <- qG*yGM/(qN*uNB*B+qG*yGM)
if(alphaVstar.1 > alphaVstar.2){ 
  alphaVstar <- alphaVstar.1
  alphaMset <- qC*(uCB*B*alphaVset2 - r)/(qG*yGM)
  linetypesetting <- 2
} else { 
    alphaVstar <- alphaVstar.2
    alphaMset <- qN*(uNB*B*alphaVset2)/(qG*yGM)
    linetypesetting <- 3
    }

alphaPstar <- 0
alphaMstar <- 1 - alphaVstar - alphaPstar


# Zero-waste lines

alphaPset <- 1-alphaMset-alphaVset2

points.to.plot <- as.data.frame(cbind(alphaPset,alphaMset,alphaVset2))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,1]>=0,]

TernaryPoints(points.to.plot[,1:3],type='l',lwd=1,lty = linetypesetting)

TernaryPoints(c(alphaPstar,alphaMstar,alphaVstar),cex=1,col='red',bg='red',pch=21)


# MIXOTROPH
TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
text(0.,0.95,'Mixotrophy Optimal',col='black')

# External resource availability
L <- 10
I <- 2.8 #3 #2
B <- 1e6

# Compute the strategy
valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){   values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }

# Add shading to the plot
# Add shading to the plot
ColourTernary(values,spectrum=colfunc(100))
#ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))

# Add analytical solution
term.a <- (qN*uNB*B-uCB*qC*B)/(qC*uCL*L-qN*uNI*I)
term.b <- (qG*yGM+uCB*qC*B)/(qG*yGM+qC*uCL*L)
term.c <- (qG*yGM+qC*r)/(qG*yGM+qC*uCL*L)
term.d <- (qC*r)/(qC*uCL*L-qN*uNI*I)

alphaVstar <- (term.c-term.d)/(term.a+term.b)
alphaPstar <- (qG*yGM*(1-alphaVstar)+qC*r-qC*uCB*B*alphaVstar)/(qG*yGM+qC*uCL*L)
alphaMstar <- 1 - alphaVstar - alphaPstar


# Zero-waste lines
alphaVset2 <- seq(from = 0, to = 1, length.out = 100)

alphaPset.1 <- (qN*uNB*B*alphaVset2-qC*(uCB*B*alphaVset2-r))/(qC*uCL*L-qN*uNI*I)
alphaPset.2 <- (qG*yGM*(1-alphaVset2)+qC*r-qC*uCB*B*alphaVset2)/(qG*yGM+qC*uCL*L)
alphaPset.3 <- (qG*yGM*(1-alphaVset2) - qN*uNB*B*alphaVset2)/(qN*uNI*I + qG*yGM)

points.to.plot <- as.data.frame(cbind(alphaPset.1,(1-alphaPset.1-alphaVset2),alphaVset2))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,1]>0 & points.to.plot[,2]>alphaMstar,]

TernaryPoints(points.to.plot[,1:3],type='l',lwd=1)

points.to.plot <- as.data.frame(cbind(alphaPset.2,(1-alphaPset.2-alphaVset2),alphaVset2))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,1]>0 & points.to.plot[,2]>0&points.to.plot[,2]<alphaMstar,]

TernaryPoints(points.to.plot[,1:3],type='l',lwd=1,lty=2)

points.to.plot <- as.data.frame(cbind(alphaPset.3,(1-alphaPset.3-alphaVset2),alphaVset2))
points.to.plot$sum <- rowSums(points.to.plot)
points.to.plot <- points.to.plot[points.to.plot[,4] <= 1 & points.to.plot[,1]>0 & points.to.plot[,2]>0&points.to.plot[,2]<alphaMstar,]

TernaryPoints(points.to.plot[,1:3],type='l',lwd=1,lty=3)

TernaryPoints(c(alphaPstar,alphaMstar,alphaVstar),cex=1,col='red',bg='red',pch=21)

Chemostat-MOCHA

First, let’s make sure we have a robust algorithm for choosing the optimal investment strategy. We’ll confirm using our ternary plot friend.

# Plot parameters
par(mar=rep(0.2, 4),mfrow=c(1,1))
colfunc <- colorRampPalette(c("black",'gray65',"white"))
resolutionsetting <- 100L


# MIXOTROPH PARAMETERS
# Respiratory costs
r <- 0.01

# Uptakes
yGM <- 1 # yield from growth
uCL <- 0.01 # carbon from light
uNI <- 1 # nutrients from inorganic nutrients
uCB <- 0.000003 # carbon from bacteria
uNB <- uCB/10 # nitrogen from bacteria

# Stoichiometries
qC <- 1
qN <- 1
qG <- 1


# External resource availability
L <- 1000
I <- 0
B <- 0.88e6 #1e6

L <- 10
I <- 2.8 #3 #2
B <- 1e6


TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))
#text(0.,0.95,'Phagotrophy Optimal',col='black')

# Compute the strategy
# Growth rate functions
jCfcn <- function (a, b, c) { qC*(uCL*L*a + uCB*B*c - r) }
jNfcn <- function(a,b,c){   qN*(uNI*I*a + uNB*B*c) }
jGfcn <- function(a,b,c){   qG*yGM*b }

valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){   values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }

# Add shading to the plot
#ColourTernary(values,spectrum=colfunc(100))
ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))

# Add analytical solution
strats <- matrix(data=NA,nrow=3,ncol=4)

## Strict phagotrophy
alpha_til_V <- max((qG*yGM+qC*r)/(qC*uCB*B + qG*yGM),(qG*yGM)/(qN*uNB*B+qG*yGM))
alpha_til_P <- 0
alpha_til_M <- 1 - alpha_til_V - alpha_til_P
g_til_CV <- qC*(uCB*B*alpha_til_V-r)
g_til_NV <- qN*uNB*B*alpha_til_V
g_til = qG*yGM*alpha_til_M # growth rate of strict phagotroph
strats[1,1] <- alpha_til_V
strats[1,2] <- alpha_til_P
strats[1,3] <- alpha_til_M
strats[1,4] <- g_til

## Strict phototrophy
alpha_bar_V = 0
alpha_bar_P = max(c((qG*yGM+qC*r)/(qC*uCL*L+qG*yGM),qG*yGM/(qN*uNI*I+qG*yGM)))
alpha_bar_M = 1-alpha_bar_P
g_bar_CP = qC*(uCL*L*alpha_bar_P - r)
g_bar_NP = qN*uNI*I*alpha_bar_P
g_bar = qG*yGM*alpha_bar_M # growth rate of strict phototroph
strats[2,1] <- alpha_bar_V
strats[2,2] <- alpha_bar_P
strats[2,3] <- alpha_bar_M
strats[2,4] <- g_bar

## Mixotrophy
term.a <- (qN*uNB*B-uCB*qC*B)/(qC*uCL*L-qN*uNI*I)
term.b <- (qG*yGM+uCB*qC*B)/(qG*yGM+qC*uCL*L)
term.c <- (qG*yGM+qC*r)/(qG*yGM+qC*uCL*L)
term.d <- (qC*r)/(qC*uCL*L-qN*uNI*I)

alpha_star_V <- (term.c-term.d)/(term.a+term.b)
alpha_star_P <- (qG*yGM*(1-alpha_star_V)+qC*r-qC*uCB*B*alpha_star_V)/(qG*yGM+qC*uCL*L)
alpha_star_M <- 1 - alpha_star_V - alpha_star_P
g_star = qG*yGM*alpha_star_M
strats[3,1] <- alpha_star_V
strats[3,2] <- alpha_star_P
strats[3,3] <- alpha_star_M
strats[3,4] <- g_star

## Picking optimum
strats_simplify <- strats[strats[,2]>=0&strats[,1]>=0,]
strats_opt <- strats_simplify[strats_simplify[,4]==max(strats_simplify[,4]),]

TernaryPoints(c(strats_opt[2],strats_opt[3],strats_opt[1]),cex=1,col='red',bg='red',pch=21)

#strats_opt

Great, so now we have an optimization. Let’s put it into a function.

mixo_opt <- function(L,I,B,qG,qC,qN,yGM,r,uCB,uNB,uCL,uNI){
  # Add analytical solution
strats <- matrix(data=NA,nrow=3,ncol=4)

## Strict phagotrophy
alpha_til_V <- max((qG*yGM+qC*r)/(qC*uCB*B + qG*yGM),(qG*yGM)/(qN*uNB*B+qG*yGM))
alpha_til_P <- 0
alpha_til_M <- 1 - alpha_til_V - alpha_til_P
g_til_CV <- qC*(uCB*B*alpha_til_V-r)
g_til_NV <- qN*uNB*B*alpha_til_V
g_til = qG*yGM*alpha_til_M # growth rate of strict phagotroph
strats[1,1] <- alpha_til_V
strats[1,2] <- alpha_til_P
strats[1,3] <- alpha_til_M
strats[1,4] <- g_til

## Strict phototrophy
alpha_bar_V = 0
alpha_bar_P = max(c((qG*yGM+qC*r)/(qC*uCL*L+qG*yGM),qG*yGM/(qN*uNI*I+qG*yGM)))
alpha_bar_M = 1-alpha_bar_P
g_bar_CP = qC*(uCL*L*alpha_bar_P - r)
g_bar_NP = qN*uNI*I*alpha_bar_P
g_bar = qG*yGM*alpha_bar_M # growth rate of strict phototroph
strats[2,1] <- alpha_bar_V
strats[2,2] <- alpha_bar_P
strats[2,3] <- alpha_bar_M
strats[2,4] <- g_bar

## Mixotrophy
term.a <- (qN*uNB*B-uCB*qC*B)/(qC*uCL*L-qN*uNI*I)
term.b <- (qG*yGM+uCB*qC*B)/(qG*yGM+qC*uCL*L)
term.c <- (qG*yGM+qC*r)/(qG*yGM+qC*uCL*L)
term.d <- (qC*r)/(qC*uCL*L-qN*uNI*I)

alpha_star_V <- (term.c-term.d)/(term.a+term.b)
alpha_star_P <- (qG*yGM*(1-alpha_star_V)+qC*r-qC*uCB*B*alpha_star_V)/(qG*yGM+qC*uCL*L)
alpha_star_M <- 1 - alpha_star_V - alpha_star_P
g_star = qG*yGM*alpha_star_M
strats[3,1] <- alpha_star_V
strats[3,2] <- alpha_star_P
strats[3,3] <- alpha_star_M
strats[3,4] <- g_star

## Picking optimum
strats_simplify <- strats[strats[,2]>=0&strats[,1]>=0,]
strats_opt <- strats_simplify[strats_simplify[,4]==max(strats_simplify[,4]),]

Res <- list("Optimum" = strats_opt)
}

# Plot parameters
par(mar=rep(0.2, 4),mfrow=c(1,1))
colfunc <- colorRampPalette(c("black",'gray65',"white"))
resolutionsetting <- 100L


# MIXOTROPH PARAMETERS
# Respiratory costs
r <- 0.01

# Uptakes
yGM <- 1 # yield from growth
uCL <- 0.01 # carbon from light
uNI <- 1 # nutrients from inorganic nutrients
uCB <- 0.000003 # carbon from bacteria
uNB <- uCB/10 # nitrogen from bacteria

# Stoichiometries
qC <- 1
qN <- 1
qG <- 1


# External resource availability
L <- 10
I <- 2.8 #3 #2
B <- 1e6


TernaryPlot(alab = 'α_P', blab = 'α_M', clab = 'α_V',axis.labels = seq(from=0,to=1,by=0.1),lab.offset=c(0.18,0.18,0.15))

# Compute the strategy for all investment strategies
jCfcn <- function (a, b, c) { qC*(uCL*L*a + uCB*B*c - r) }
jNfcn <- function(a,b,c){   qN*(uNI*I*a + uNB*B*c) }
jGfcn <- function(a,b,c){   qG*yGM*b }

valuesC <- TernaryPointValues(jCfcn, resolution=resolutionsetting)
valuesN <- TernaryPointValues(jNfcn, resolution=resolutionsetting)
valuesG <- TernaryPointValues(jGfcn, resolution=resolutionsetting)
values <- valuesC
for(i in 1:dim(valuesC)[2]){   values[3,i] <- min(valuesC[3,i],valuesN[3,i],valuesG[3,i]) }

# Add shading to the plot
#ColourTernary(values,spectrum=colfunc(100))
ColourTernary(values,spectrum=viridisLite::viridis(500L,alpha=.8))

# Compute the growth-maximizing strategy
Res <- mixo_opt(L,I,B,qG,qC,qN,yGM,r,uCB,uNB,uCL,uNI)

TernaryPoints(c(Res$Optimum[2],Res$Optimum[3],Res$Optimum[1]),cex=1,col='red',bg='red',pch=21)

Without dynamic resources, if the dilution rate is lower than the initial growth rate, the population will grow without bound.

# Parameters
# External resource availability
L0 <- 100
I0 <- 10
B0 <- 1e6

# Respiratory costs
r <- 0.1#0.01

# Uptakes
yGM <- 1
uCL <- 0.01
uNI <- 1
uCB <- 0.000003
uNB <- uCB/10

# Stoichiometries
qC <- 1
qN <- 1
qG <- 1

# Dilution rate of the chemostat
D <- 0.001

# Set up the simulation
tset <- seq(from = 0, to = 1000, by = 0.01)
Res <- mixo_opt(L,I,B,qG,qC,qN,yGM,r,uCB,uNB,uCL,uNI)
Vset <- NaN*tset; Vset[1] <- Res$Optimum[1]
Pset <- Vset; Pset[1] <- Res$Optimum[2]
Mset <- Vset; Mset[1] <- Res$Optimum[3]
Bset <- Vset; Bset[1] <- B0
Iset <- Vset; Iset[1] <- I0
Lset <- Vset; Lset[1] <- L0

growth <- NaN*tset; growth[1] <- Res$Optimum[4]

# Run the simulation
for(i in 2:length(tset)){
  dt <- tset[i] - tset[i-1]
  
  L <- L0
  B <- B0
  I <- I0
  
  jC <- ( uCL*L*Pset[i-1] + uCB*B*Vset[i-1] - r*(Pset[i-1]+Vset[i-1]+Mset[i-1]) )
  jN <- ( uNI*I*Pset[i-1] + uNB*B*Vset[i-1] )
  jG <- ( yGM*Mset[i-1] )
  
  growth[i] <- min(qC*jC,qN*jN,qG*jG)
  
  Res <- mixo_opt(L,I,B,qG,qC,qN,yGM,r,uCB,uNB,uCL,uNI)
  alphaV <- Res$Optimum[1]
  alphaP <- Res$Optimum[2]
  alphaM <- Res$Optimum[3]
  
  Pset[i] <- Pset[i-1] + dt*( growth[i]*alphaP - D*Pset[i-1])
  Vset[i] <- Vset[i-1] + dt*( growth[i]*alphaV - D*Vset[i-1])
  Mset[i] <- Mset[i-1] + dt*( growth[i]*alphaM - D*Mset[i-1])
    
}


Tset <- Pset+Vset+Mset

plot(tset,Pset/Tset,col=Pcol,lwd=2,ylim=c(0,1),type='l',las=1,xlab='',ylab='')
abline(h=alphaP,lty=2,col=Pcol)
lines(tset,Vset/Tset,col=Vcol,lwd=2)
abline(h=alphaV,col=Vcol,lty=2)
lines(tset,Mset/Tset,col=Mcol,lwd=2)
abline(h=alphaM,col=Mcol,lty=2)
mtext('Time',side=1,line=2,cex=0.9)
mtext('Proportion of Biomass',side=2,line=2.5,cex=0.9)
legend(x='top',inset=0.05,legend=c('Vacuoles','Plastids','Machinery'),lwd=2,col=c(Vcol,Pcol,Mcol),horiz=T,cex=.8)
mtext(c('α_V','α_P','α_M'),side=4, line = 0.2, at = c(alphaV,alphaP,alphaM),las=1,cex=0.8)

plot(tset,Pset,col=Pcol,lwd=2,type='l',las=1,xlab='',ylab='')
#abline(h=alphaP,lty=2,col=Pcol)
lines(tset,Vset,col=Vcol,lwd=2)
#abline(h=alphaV,col=Vcol,lty=2)
lines(tset,Mset,col=Mcol,lwd=2)
#abline(h=alphaM,col=Mcol,lty=2)
mtext('Time',side=1,line=2,cex=0.9)
mtext('Biomass',side=2,line=2.5,cex=0.9)
legend(x='top',inset=0.05,legend=c('Vacuoles','Plastids','Machinery'),lwd=2,col=c(Vcol,Pcol,Mcol),horiz=T,cex=.8)

#mtext(c('α_V','α_P','α_M'),side=4, line = 0.2, at = c(alphaV,alphaP,alphaM),las=1,cex=0.8)

Adding in chemostat dynamics for the mixotroph

# Parameters
# External resource availability
L0 <- 100
I0 <- 10
B0 <- 1e6

# Respiratory costs
r <- 0.1#0.01

# Uptakes
yGM <- 1
uCL <- 0.01
uNI <- 1
uCB <- 0.000003
uNB <- uCB/10

# Stoichiometries
qC <- 1
qN <- 1
qG <- 1

# Dilution rate of the chemostat
D <- 0.001

# Set up the simulation
tset <- seq(from = 0, to = 1000, by = 0.01)
Res <- mixo_opt(L0,I0,B0,qG,qC,qN,yGM,r,uCB,uNB,uCL,uNI)
Vset <- NaN*tset; Vset[1] <- Res$Optimum[1]*.5
Pset <- Vset; Pset[1] <- Res$Optimum[2]
Mset <- Vset; Mset[1] <- Res$Optimum[3]
Bset <- Vset; Bset[1] <- B0
Iset <- Vset; Iset[1] <- I0
Lset <- Vset; Lset[1] <- L0

growth <- NaN*tset; growth[1] <- Res$Optimum[4]

# Run the simulation
for(i in 2:length(tset)){
  dt <- tset[i] - tset[i-1]
  
  L <- L0
  B <- B0
  I <- I0
  
  jC <- ( uCL*L*Pset[i-1] + uCB*B*Vset[i-1] - r*(Pset[i-1]+Vset[i-1]+Mset[i-1]) )
  jN <- ( uNI*I*Pset[i-1] + uNB*B*Vset[i-1] )
  jG <- ( yGM*Mset[i-1] )
  
  growth[i] <- min(qC*jC,qN*jN,qG*jG)
  
  Res <- mixo_opt(L,I,B,qG,qC,qN,yGM,r,uCB,uNB,uCL,uNI)
  alphaV <- Res$Optimum[1]
  alphaP <- Res$Optimum[2]
  alphaM <- Res$Optimum[3]
  
  Pset[i] <- Pset[i-1] + dt*( growth[i]*alphaP - D*Pset[i-1])
  Vset[i] <- Vset[i-1] + dt*( growth[i]*alphaV - D*Vset[i-1])
  Mset[i] <- Mset[i-1] + dt*( growth[i]*alphaM - D*Mset[i-1])
    
}


Tset <- Pset+Vset+Mset

par(mar=c(1,4,3,2),mfrow=c(2,1))

plot(tset,Pset/Tset,col=Pcol,lwd=2,ylim=c(0,1),type='l',las=1,xlab='',ylab='',xaxt='n')
axis(side=1,labels=F)
abline(h=alphaP,lty=2,col=Pcol)
lines(tset,Vset/Tset,col=Vcol,lwd=2)
abline(h=alphaV,col=Vcol,lty=2)
lines(tset,Mset/Tset,col=Mcol,lwd=2)
abline(h=alphaM,col=Mcol,lty=2)
#mtext('Time',side=1,line=2,cex=0.9)
mtext('Proportion of Biomass',side=2,line=2.5,cex=0.9)
legend(x='top',inset=0.05,legend=c('Vacuoles','Plastids','Machinery'),lwd=2,col=c(Vcol,Pcol,Mcol),horiz=T,cex=.8)
mtext(c('α_V','α_P','α_M'),side=4, line = 0.2, at = c(alphaV,alphaP,alphaM),las=1,cex=0.8)

par(mar=c(4,4,0,2))
plot(tset,Pset,col=Pcol,lwd=2,type='l',las=1,xlab='',ylab='')
#abline(h=alphaP,lty=2,col=Pcol)
lines(tset,Vset,col=Vcol,lwd=2)
#abline(h=alphaV,col=Vcol,lty=2)
lines(tset,Mset,col=Mcol,lwd=2)
#abline(h=alphaM,col=Mcol,lty=2)
mtext('Time',side=1,line=2,cex=0.9)
mtext('Biomass',side=2,line=2.5,cex=0.9)
legend(x='top',inset=0.05,legend=c('Vacuoles','Plastids','Machinery'),lwd=2,col=c(Vcol,Pcol,Mcol),horiz=T,cex=.8)

#mtext(c('α_V','α_P','α_M'),side=4, line = 0.2, at = c(alphaV,alphaP,alphaM),las=1,cex=0.8)

And now dynamic resources

# Parameters
# External resource availability
L0 <- 100 #1000
I0 <- 10 #0#10
B0 <- 1e6

# Respiratory costs
r <- 1#0.1#0.01

# Uptakes
yGM <- 1
uCL <- 0.01 #0.01*2
uNI <- 1
attack_rate <- 0.000003
f <- 1 # conversion efficiency of bacterial C to mixotroph C
uCB <- attack_rate*f
uNB <- uCB/10 ##5*attack_rate*f

# Stoichiometries
qC <- 1
qN <- 1
qG <- 1

# Dilution rate of the chemostat
D <- 0.1

# Absorptivities of chemostat components
kappa0 <- 0 # base absorptivity of water
kappaI <- 0 # absorptivity of inorganic nutrients
kappaB <- 0.#00001 # per capita absorptivity of bacteria
kappaP <- 0.#001 # per capita absorptivity of plastids
kappaV <- 0.#00001 # per capita absorptivity of vacuoles
kappaM <- 0.#00001 # per capita absorptivity of growth machinery

# Set up the simulation
tset <- seq(from = 0, to = 20, by = 0.001)
Res <- mixo_opt(L0,I0,B0,qG,qC,qN,yGM,r,uCB,uNB,uCL,uNI)
Vset <- NaN*tset; Vset[1] <- Res$Optimum[1]
Pset <- Vset; Pset[1] <- Res$Optimum[2]
Mset <- Vset; Mset[1] <- Res$Optimum[3]
Bset <- Vset; Bset[1] <- B0
Iset <- Vset; Iset[1] <- I0
Lset <- Vset; Lset[1] <- L0
alphaVset <- Vset; alphaVset[1] <- Res$Optimum[1]
alphaPset <- Vset; alphaPset[1] <- Res$Optimum[2]
alphaMset <- Vset; alphaMset[1] <- Res$Optimum[3]
growth <- NaN*tset; growth[1] <- Res$Optimum[4]

# Run the simulation
for(i in 2:length(tset)){
  dt <- tset[i] - tset[i-1]
  
  # Determine resource acquisition
  L <- Lset[i-1]
  B <- Bset[i-1]
  I <- Iset[i-1]
  jC <- ( uCL*L*Pset[i-1] + uCB*B*Vset[i-1] - r*(Pset[i-1]+Vset[i-1]+Mset[i-1]) )
  jN <- ( uNI*I*Pset[i-1] + uNB*B*Vset[i-1] )
  jG <- ( yGM*Mset[i-1] )
  
  # Compute growth rate
  growth[i] <- min(qC*jC,qN*jN,qG*jG)
  
  # Determine optimal investment strategy
  Res <- mixo_opt(L,I,B,qG,qC,qN,yGM,r,uCB,uNB,uCL,uNI)
  alphaVset[i] <- Res$Optimum[1]
  alphaPset[i] <- Res$Optimum[2]
  alphaMset[i] <- Res$Optimum[3]
  
  # Update mixotroph state variables
  Pset[i] <- Pset[i-1] + dt*( growth[i]*alphaPset[i] - D*Pset[i-1])
  Vset[i] <- Vset[i-1] + dt*( growth[i]*alphaVset[i] - D*Vset[i-1])
  Mset[i] <- Mset[i-1] + dt*( growth[i]*alphaMset[i] - D*Mset[i-1])
  
  # Update resource state variables
  Bset[i] <- Bset[i-1] + dt*( D*(B0 - Bset[i-1]) - attack_rate*Vset[i-1]*Bset[i-1] )
  Iset[i] <- Iset[i-1] + dt*( D*(I0 - Iset[i-1]) - uNI*Iset[i-1]*Pset[i-1] )
  Lset[i] <- L0*exp(-kappa0 - kappaP*Pset[i] - kappaV*Vset[i] - kappaM*Mset[i] - kappaB*Bset[i] - kappaI*Iset[i])
    
}


Tset <- Pset+Vset+Mset

par(mar=c(1,4,3,2),mfrow=c(2,1))

plot(tset,Pset/Tset,col=Pcol,lwd=2,ylim=c(0,1),type='l',las=1,xlab='',ylab='',xaxt='n')
axis(side=1,labels=F)
#abline(h=alphaP,lty=2,col=Pcol)
lines(tset, alphaPset, lty = 2, col = Pcol)
lines(tset,Vset/Tset,col=Vcol,lwd=2)
#abline(h=alphaV,col=Vcol,lty=2)
lines(tset,alphaVset, lty = 2, col = Vcol)
lines(tset,Mset/Tset,col=Mcol,lwd=2)
#abline(h=alphaM,col=Mcol,lty=2)
lines(tset,alphaMset,lty=2, col = Mcol)
#mtext('Time',side=1,line=2,cex=0.9)
mtext('Proportion of Biomass',side=2,line=2.5,cex=0.9)
legend(x='right',inset=0.05,legend=c('Vacuoles','Plastids','Machinery'),lwd=2,col=c(Vcol,Pcol,Mcol),horiz=F,cex=.8)
mtext(c('α_V','α_P','α_M'),side=4, line = 0.2, at = c(alphaV,alphaP,alphaM),las=1,cex=0.8)

par(mar=c(4,4,0,2))
plot(tset,Pset,col=Pcol,lwd=2,type='l',las=1,xlab='',ylab='',ylim=c(min(Pset,Vset,Mset,Tset),max(Pset,Vset,Mset,Tset)),log='y')
#abline(h=alphaP,lty=2,col=Pcol)
lines(tset,Vset,col=Vcol,lwd=2)
#abline(h=alphaV,col=Vcol,lty=2)
lines(tset,Mset,col=Mcol,lwd=2)
#abline(h=alphaM,col=Mcol,lty=2)
mtext('Time',side=1,line=2,cex=0.9)
mtext('Biomass',side=2,line=2.5,cex=0.9)
lines(tset,Tset,col='black',lwd=2)

#legend(x='top',inset=0.05,legend=c('Vacuoles','Plastids','Machinery'),lwd=2,col=c(Vcol,Pcol,Mcol),horiz=T,cex=.8)
#mtext(c('α_V','α_P','α_M'),side=4, line = 0.2, at = c(alphaV,alphaP,alphaM),las=1,cex=0.8)

plot(tset,alphaPset,type='l')

plot(Bset,Iset,type='l',lwd=2)